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Abstract. - A lattice-Boltzmann model for the study of the dynamics of oil-water-surfactant 
mixtures is constructed. The model, which is based on a Ginzburg-Landau theory of amphiphilic 
systems with a single, scalar order parameter, is then used to calculate the spectrum of undula- 
tion modes of an oil-water interface and the spontaneous emulsification of oil and water after a 
quench from two-phase coexistence into the lamellar phase. A comparison with some analytical 
results shows that the model provides an accurate description of the static and dynamic behavior 
of amphiphilic systems. 



Introduction. - Self-assembling amphiphilic system are characterized by a large variety of 
phases jl], ||. Many of these are structured on mesoscopic length scales. Examples include 
the swollen lamellar phase, the microcmulsion phase in ternary amphiphilic mixtures, and 
the sponge phase in aqueous amphiphile solutions. In all of these cases, the amphiphile 
assembles into two-dimensional sheets which separate space into two networks of oil- and/or 
water-channels, both of which span the whole system. 

The mesoscopic structure leads to dynamical behavior which differs qualitatively from that 
of simple fluid mixtures. In addition to the diffusive dynamics of all components, hydrodynamic 
flow plays an important role. Three classes of theoretical models which include hydrodynamic 
effects have been employed so far. The first are time-dependent Ginzburg-Landau models 
[PI El' D w hich are coupled to a linearized Navier-Stokes equation (i.e. generalizations of model 
H of Hohenberg and Halperin [|| to amphiphilic systems). The second are membrane models, in 
which the hydrodynamics of the solvent is taken into account || . Finally, a lattice-gas model 

[l^] — similar in spirit to vector lattice models employed to study equilibrium properties of 
microemulsions jl| — has been formulated and applied to study the dynamics of these systems. 

An interesting alternative to these descriptions of the hydrodynamics of complex fluids is 
the Lattice-Boltzmann (LB) approach [O, O, 13, It has recently been shown how LB 
methods can be generalized to study the dynamics of non-ideal fluids p"5| , |l6| ] and fluid mixtures 
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[ pL T[ , [L6|, [L8| with a given free-energy functional. In this paper we show how this approach can 
be extended to study the dynamical behavior of amphiphilic systems. 

Lattice- Boltzmann Approach. - The basic variables of a LB model for two-component fluid 
mixtures are the "microscopic" distribution functions fi(r) and g'i(r), which are related to the 
density p, density difference Ap, and velocity u by 

p(r)=J2fi(r); Ap(r) = £>(r); p(r)u(r) = £ ftr)* , (1) 

i i i 

where i denotes the possible velocity states, eo = 0, and with i > 1 are the unit lattice 
vectors. The dynamics of the system is determined by the temporal evolution equations 



/.(r + e^f+l) = f i (r,t)--[f i (r,t)-f^(r,t)] (2) 
ft (r + ei,i+l) - g i {v,t)-—\ gi {v,t)-gl q {v,i)}. (3) 

TA 

It has been shown in rcf. |l7| that the time dependence of Ap(r), p(r), and the momentum 
density is described by Cahn-Hilliard and Navier-Stokes equations on sufficiently large length 
scales when the zeroth and first moments of f^ 9 and g\ q are required to fulfill eq. (Q), and, in 
addition 

^2 fi" e iai e i0 = P a(3 + PU a Up (4) 

i 

2j#r e ia e i/3 = r Ap,6 a p + Apu a Up , (5) 

i 

where A/i = 5J-[p, Ap]/5Ap is the chemical potential for the density difference Ap, T is the 
mobility, and P a p the pressure tensor. We choose ta = (1 + 1/ v3) /2 to minimize terms which 
break Galilean invariance |l7|| . 

Amphiphilic Systems. - Many equilibrium properties of ternary amphiphilic systems can 
be accurately described using a Ginzburg-Landau model for a single, scalar order parameter 
<fi(r), which we identified with the local density difference, Ap(r), of oil and water. This model 
is defined by the free energy functional [|l9| |J 

= J d d r{c(V 2 0) 2 + 3 (^)(V0) 2 + /(</>)} . (6) 

We assume that the full free-energy functional, which includes the density field, can be written 
as T[p,4>\ — J-<j>[4>] + Fplp], where T p \p\ — J d d r{biy p) 2 4- ip(p)}. A similar model, with 
5(0) = .9o = const, in eq. (||) and a different ansatz for the free energy density f(cj>) (see 
below), has been employed in ref. |Q to study the spinodal decomposition to the lamellar 
phase. 

The first task in the construction of the LB model is to determine the pressure tensor P a /3, 
which is defined by the condition of hydrostatic equilibrium 

d P af} (r) = p(r)d a p{r) + Ap(r)d a Ap{r) , (7) 

in a system with external potential /i(r) for the density p(r). The pressure tensor is not unique, 
since eq. (Q) is invariant under the substitution P a p — > P a p + G a p with dpG a p = 0. Indeed, a 
Chapman-Enskog expansion shows that only the divergence of the pressure tensor appears in 
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the Cahn-Hilliard and Navier-Stokes equations. Eq. (Q) can be solved by expanding in powers 
of 4>(r) = Ap(r) and its derivatives. After some algebra, we obtain 

Paf3 = PoSa/3 + Pa/3 , (8) 

where 

po = 2 C 0(V 2 V^)-c(V 2 0) 2 -2 5 (0) ( /)V 2 0- 5 (0)(V0) 2 -.g'(0) ( /)(V ( /)) 2 

-2bpV 2 p - b(Vp) 2 + #'(0) - f(<(>) + (nP'(p) - r/>{p) (9) 

is the (isotropic) local thermodynamic pressure, and 

P a0 = 2b(d aP ) (d l3 p) + 2g(4,)(d a (f>) {d p 4>) +cK{V 2 4>) 2 5 aP -2c{d a d^){d fi d^) 

+ C (A - 2){{d a 4>)(d p v 2 <p) + {di3<p){d a v 2 <p)} + 2c(i - AXSafya^xa^) 

+cA(f > {d a d l 3V 2 (f>) - cK4>{V 2 V 2 4>)8 a p + c(4 - A)(5 ct 5 (3 0)(V 2 0) , (10) 

with an arbitrary real number A, which reflects the non-uniqueness of the pressure tensor. 
For c = and g((f>) = go eq. ( |Io|) reduces to the well known expression for inhomogeneous 
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Since this result is rather lengthy, it is important to see that eqs. (|l|) and (^) can be solved 
explicitly for any pressure tensor. In the case of the two-dimensional triangular lattice, the 
equilibrium distribution f eq is found to be 

]l q = P -TrP afj -pu 2 (11) 



/f - - \TrP a p + | (e, • u) - ^u 2 + |p(ej • u) 2 + ^P af3 e ia e i0 . (12) 



The distribution g eq can be obtained from eqs. ( |11|) and ( J12[ ) for f eq by replacing P a p by 
T(6J-[p, Ap]/6 Ap) Sap. It can be shown easily that this result agrees with the solutions of the 
special cases studied in refs. ]l5|, |l7], 

Equilibrium Properties. - To test the correctness of the approach described above, we 



have followed refs. fllq, 17, [16| and first studied the equilibrium properties of our model. We 



have checked both the density profiles and the phase diagram for the Ginzburg-Landau model 
introduced in ref. |21j] , with f(4>) — cj (</ )2 — ^ 2 ) 2 (</ )2 + fo) and g((f>) = g + 32</> 2 , and for a 
simple ideal-gas free-energy density ip(p) = Tp\n(p). The deviation of the LB results from 
the solutions of the mean-field equations is on the order of a few percent if the parameters 
are chosen such that the interface width is at least several lattice constants. The phase 
diagram is determined by studying the stability of the microemulsion phase with respect to 
small density fluctuations. Similarly, the phase transition between the lamellar phase and 
oil-water coexistence is located by looking for either the growth of lamellae or a coarsening 
of the bicontinuous oil-water network after a quench from the microemulsion phase. In both 
cases, the transitions are found to agree within a few percent with those of the mean-field 
phase diagram. 

Undulation Modes. - Due to the presence of the amphiphile, interfaces in amphiphilic 
systems can have an ultralow interfacial tension a. The fluctuations of an interface are 
controlled in this limit by the bending rigidity k. For the typical length scales of amphiphilic 
systems, the hydrodynamic flow is usually in the low- Reynolds-number, creeping-flow regime. 
In this case, a calculation of membrane fluctuations in the Stokes approximation [Q shows 
that the spectrum of the undulation modes of wavenumber k is given by 

w(jfe) = ±(a + Kk 2 ^jk , (13) 
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Fig. 1. - Spectrum of undulation modes. The parameters are c = 0.01, u>b = 0.01, (j>b = 0.5, fo = 0.1, 
g 2 = 4(j} b {y/cuj(l + 4/ ) - g ) + 0.01, A = 0.123, b = 0.05, T = 0.5, V = 0.1, and r p = 5.0. Results are 
shown for go = +0.004 (o, dashed line) and go = —0.004 (o, full line), which correspond to large and 
small a, respectively. 

Fig. 2. - Interface width w after a quench of a system from oil-water coexistence into the lamellar 
phase. The parameters are the same as in fig. [j], except for t p = 1.0 and go — —0.0125. The symbols 
indicate times for which the interfacial tension profile is shown in fig. 0. 



where r\ is the shear viscosity. We have determined the relaxation times of the undulation 
modes of different wavelengths from our LB model of lattice sizes 64 x 128 and 128 x 128 with 
periodic boundary conditions. The results are shown in fig. [l] for two values of the parameter 
go, which we use to tune the interfacial tension One of the big advantages of the LB 

approach used here is that all the coefficients in eq. ( |i~3| ) can be determined independently 
- a and n from a study of the elasticity of curved interfaces in the Ginzburg-Landau model 
and 77 = (2r p -l)/8fr om a Chapman-Enskog expansion of the Boltzmann equation (f3) 
[16 . The simulation results can therefore be compared directly with theory. It can be seen 
from fig. |l| for two values of the interfacial tension that the agreement is very good; the almost 
perfect agreement for go = —0.004 is fortuitous, however. 

We have checked explicitly that our results do not depend on the parameter A in the pressure 
tensor, eq. (|i"o|), as they should. 

Spontaneous Emulsification. - When a system is taken abruptly from oil-water coexistence 
into the region of the phase diagram where the lamellar phase is stable, the large oil- and 
water-domains have to break up and form a intertwined structure on a much smaller length 
scale. To study this process in the LB scheme, we equilibrate a system of equal oil- and 
water-content in a metastable state of oil- water coexistence — with a planar interface — inside 
the lamellar regime. We then disturb the interface position by a sinusoidal capillary wave of the 
largest possible wavelength. The behavior of the interface is monitored by measuring its width 
w = (1 /L) -\/2^j h? (xi), where L is the number of lattice points, and x is the coordinate parallel 
to, and h(x) is the local distance from, the initial planar interface. The time dependence of 
w is shown in fig. |^. The width first decreases due to a relaxation of the local interface 
profile, and a decrease in the amplitude of the capillary wave. After about 10 5 time steps, 
w increases again due to the increase of the amplitude of a mode of half the wavelength 
of the original perturbation. This mode is excited due to anharmonic terms in an effective 
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Fig. 3. - Typical flow pattern some time after the quench, where a large deformation of the interface 
has developed. Only a part of the full lattice and a subset of the velocity vectors is shown. The 
parameters are the same as in fig. El 



Fig. 4. - (a) Variation of the local interfacial tension a(x 
after the quench (from top to bottom), as indicated in fig^ 
time shown in (a) . The parameters are the same as in fig. |2 



parallel to the interface at various times 
L (b) Interface position h(x) at the last 



interface Hamiltonian which lead to the coupling of different — in particular to neighboring 
undulation modes. The new mode dominates because its characteristic frequency is larger 
than that of the original mode, compare eq. This mode grows slowly until, after about 

6 • 10 5 time steps, w becomes comparable to the interface width. A typical flow pattern at 
this time is shown in fig. |[ At later stages of emulsification, the interface shape becomes very 
complicated, and a characterization by the width w is no longer possible. 

It is interesting to study the early stages of the emulsification process in more detail. To 
do so, we have calculated the local interfacial tension, which is given by 

a(x) = J dz[P N (v)-P T (r)] , (14) 

where Pjv = P zz and Pt = P xx are the normal and tangential components of the pressure 
tensor, respectively. The result is shown in fig. |]for different times after the quench. After the 
initial relaxation, the interfacial tension takes a value, which is positive, but nearly vanishes. 
After 2 • 10 5 time steps, a periodic structure appears, with a wavelength half of that of the 
undulation mode. The smallest tensions are found at the maxima of the function h(x), i.e. 
where the interface motion is the strongest, while the largest tensions are located at the 
inflection points of h(x). After about 4 • 10 5 time steps, o~(x) becomes negative at its minima. 
The growth of the interface width now accelerates; this is accompanied by a further decrease 
of the interfacial tension. 

Thus, we find that the surface tension along the interface becomes inhomogeneous in the 
early stages of the emulsification process. This indicates that most of the new interfacial area 
is created at the points where the interface moves fastest. 

Summary and Conclusions. We have shown in this paper how the LB scheme for 
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non-ideal fluid mixtures, which has been introduced in refs. |17| , [16| , can be extended to 
describe amphophilic systems. We have tested both the equilibrium and dynamical properties 
of our model and have found excellent agreement of the results with complementary theoretical 
approaches. The model has then been used to study the spontaneous emulsification of oil and 
water after a quench into the lamellar phase. 

Another conserved density, which we haven't taken into account in our present model, is the 
density of the amphiphile itself. Ginzburg-Landau models, which include this density, have 
already been investigated in thermal equilibrium Q , so that the construction of a LB model 
for this case should be straightforward. The diffusion of the amphiphile to the interface should 
lead to new dynamical phenomena [ p2[ . 
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